!     ---------------
!     PROGRAM IPS2000
!     ---------------
!
!     Example for using the subroutine IPS2000:
!     - Reading the ASCII planetary files.
!     - Computation of planetary rectangular coordinates.
!
!     First Julian date: JD1903670.5.
!     Step length: 73000.0 days.
!     Number of dates: 11.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*10 pla(8)
      dimension r(3)
      data pla/
     .'MERCURY','VENUS','EARTH-MOON','MARS','JUPITER',
     .'SATURN','URANUS','NEPTUNE'/
!
! --- Parameters -------------------------------------------------------
!
      t0=1903670.5d0
      step=73000.0d0
      ndat=11
      lu=10
!
! --- Results file (unit 20) -------------------------------------------
!
      open (20,file='ips2000.txt')
      write (20,1001)
      write (*,1001)
!
! --- Compute the planetary coordinates --------------------------------
!
      do ipla=1,8
         t=t0
         write (20,1002) pla(ipla)
         write (*,1002)  pla(ipla)         
         do n=1,ndat
            call IPS2000 (t,ipla,lu,r,ierr)
            if (ierr.ne.0) then
               write (*,1003) ierr
               stop
            endif
            write (20,1004) t,r
            write (*,1004)  t,r
            t=t+step
         enddo
         if (ipla.ne.8) then
            write (*,1005)
            read (*,*)
         else
            write (*,1006)
         endif
      enddo
      stop
!
! --- Formats ----------------------------------------------------------
!
1001  format (/2x,'IPS2000 : PLANETARY COORDINATES'//
     .         2x,'heliocentric rectangular coordinates',
     .         2x,'(ecliptic and equinox of J2000)')
1002  format (/2x,a/)     
1003  format (/2x,'ERREUR : ',i2/)
1004  format (2x,'JD',f9.1,
     .        2x,'X = ',f14.10,2x,'Y :',f14.10,2x,'Z :',f14.10,2x,'au')
1005  format (/2x,'Hit Enter ... ')
1006  format (/2x,' Done '/)
!
      end
!
!
!
      subroutine IPS2000 (tjd,ipla,lu,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : JCGF0201
!
! --- Object -----------------------------------------------------------
!
!     Improved Planetary Solutions over 2000 years.
!     Planetary heliocentric rectangular coordinates.
!
!     Source   : Planetary series IPS2000 (Improved Planetary Solutions)
!     Frame    : Dynamical equinox and ecliptic of J2000
!     Validity : 20/12/499 (JD1903670.5) - 23/8/2498 (JD2633670.5)
!
! --- Input ------------------------------------------------------------
!
!     tjd       Julian Date TDB (double precision).
!
!     ipla      Planet index (integer).
!               1 : Mercury                 5 : Jupiter
!               2 : Venus                   6 : Saturn
!               3 : Earth-Moon barycenter   7 : Uranus
!               4 : Mars                    8 : Neptune
!
!     lu        Logical unit for the IPS2000 files (integer).
!
! --- Output -----------------------------------------------------------
!
!     r(3)      Rectangular coordinates (double precision).
!               r(1) : X component (au).
!               r(2) : Y component (au).
!               r(3) : Z component (au).
!
!     ierr      Error index (integer).
!               1 : Date error.
!               2 : Planet error.
!               3 : File error.
!
! --- Remark -----------------------------------------------------------
!
!     The IPS2000 files (logical unit lu) are in the current directory:
!     ips2000.mer : Mercury       ips2000.jup : Jupiter
!     ips2000.ven : Venus         ips2000.sat : Saturn
!     ips2000.emb : Earth-Moon    ips2000.ura : Uranus
!     ips2000.mar : Mars          ips2000.nep : Neptune
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*3 ext(8)
      character*40 fname
      logical first
!
      dimension r(3)
      dimension t(0:5)
      dimension a(40000),b(40000),c(40000)
      dimension maxo(8,3),maxt(8,3,0:5),mint(8,3,0:5)
!
      save
!      
      data first/.true./
!
      data t2000/2451545.0d0/
!
      data tmin/1903670.5d0/
      data tmax/2633670.5d0/
!
      data ext/'mer','ven','emb','mar','jup','sat','ura','nep'/
!
! --- Initialization ---------------------------------------------------
!
      do i=1,3
         r(i)=0.d0
      enddo
!
      t(0)=1.d0
      t(1)=(tjd-t2000)/365250.d0
      do n=2,5
         t(n)=t(n-1)*t(1)
      enddo
!
! --- Check the parameters ---------------------------------------------
!
      if (tjd.lt.tmin.or.tjd.gt.tmax) then
         ierr=1
         return
      endif
!
      if (ipla.lt.1.or.ipla.gt.8) then
         ierr=2
         return
      endif
!
! --- Files reading ----------------------------------------------------
!
      if (first) then
!
         nr=0
!
         do ip=1,8
!
            fname='ips2000.'//ext(ip)
            open (lu,file=fname,err=200)
            read (lu,1001,err=200)
!
            do
              read (lu,1002,end=100,err=200) iv,io,nt
              maxo(ip,iv)=io
              mint(ip,iv,io)=nr+1
              maxt(ip,iv,io)=nr+nt
              do n=mint(ip,iv,io),maxt(ip,iv,io)
                 read (lu,1003,err=200) a(n),b(n),c(n)
              enddo
              nr=nr+nt
            enddo
!
100         continue
            close (lu)
!
         enddo
!
         first=.false.
!
      endif
!
! --- Compute the coordinates of the planet ----------------------------
!
      do iv=1,3
         max=maxo(ipla,iv)
         do io=0,max
         do n=mint(ipla,iv,io),maxt(ipla,iv,io)
            r(iv)=r(iv)+a(n)*t(io)*sin(b(n)*t(1)+c(n))
         enddo
         enddo
      enddo
!
      ierr=0
      return
!
! --- File error -------------------------------------------------------
!
200   continue
!
      ierr=3
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (1x)
1002  format (37x,i1,12x,i1,12x,i4)
1003  format (5x,3f25.16)
!
      return
      end
